高齢化のピークはいつか?(市区町村別の簡易推定の試行)

元データのダウンロード・整形
library(tidyverse)
library(fs)
library(readxl)
library(arrow)

dir_create("when-is-the-peak-of-the-aging-population")

file_agingrate <- "when-is-the-peak-of-the-aging-population/agingrate.xls"
if (!file_exists(file_agingrate)) {
  curl::curl_download("https://www.ipss.go.jp/pp-shicyoson/j/shicyoson18/2gaiyo_hyo/kekkahyo3_3.xls",
                      "when-is-the-peak-of-the-aging-population/agingrate.xls") 
}

agingrate <- read_excel("when-is-the-peak-of-the-aging-population/agingrate.xls",
                        range = "A4:K1861") |> 
  rename_with(\(x) c("city_code", "city_type", "pref_name", "city_name"),
              1:4) |> 
  pivot_longer(!starts_with(c("city", "pref")),
               names_to = "year",
               names_transform = list(year = parse_number),
               values_to = "agingrate") |> 
  mutate(city_code = city_code |> 
           str_pad(5,
                   pad = "0"),
         agingrate = agingrate |> 
           units::set_units(`%`) |> 
           units::set_units(1) |> 
           units::drop_units())

write_parquet(agingrate, "when-is-the-peak-of-the-aging-population/agingrate.parquet")
パッケージのロード
library(tidyverse)
library(arrow)
library(leaflet)
library(htmltools)

library(jpcity)

theme_set(theme_light())
データ読込み・傾きの計算
agingrate <- read_parquet("when-is-the-peak-of-the-aging-population/agingrate.parquet") |> 
  select(city_code, year, agingrate)

agingrate_logit <- agingrate |> 
  mutate(agingrate_logit = gtools::logit(agingrate),
         .keep = "unused")

agingrate_logit_diff <- agingrate_logit |>  
  mutate(across(c(year, agingrate_logit),
                list(lag = lag)),
         .by = city_code) |> 
  drop_na(ends_with("lag")) |> 
  mutate(year = (year + year_lag) / 2,
         agingrate_logit_diff = (agingrate_logit - agingrate_logit_lag) / (year - year_lag),
         .keep = "unused")

agingrate_logit_diff_diff <- agingrate_logit_diff |> 
  mutate(across(c(year, agingrate_logit_diff),
                list(lag = lag)),
         .by = city_code) |> 
  drop_na(ends_with("lag")) |> 
  mutate(year = (year + year_lag) / 2,
         agingrate_logit_diff_diff = (agingrate_logit_diff - agingrate_logit_diff_lag) / (year - year_lag),
         .keep = "unused")
t検定
t_test_agingrate_logit_diff_diff <- agingrate_logit_diff_diff |> 
  infer::t_test(response = agingrate_logit_diff_diff)
knitr::kable(t_test_agingrate_logit_diff_diff)
t検定の結果
statistic t_df p_value alternative estimate lower_ci upper_ci
-29.46549 9284 0 two.sided -0.0023247 -0.0024793 -0.00217
ピークの位置と値の算出
quantile_agingrate <- function(x) {
  quantile(x, c(0.025, 0.975)) |> 
    bind_rows() |> 
    rename_with(\(x) c("lower", "upper"))
}

slope_agingrate <- t_test_agingrate_logit_diff_diff$estimate / 2
axis_agingrate <- agingrate_logit_diff |> 
  mutate(axis = year - agingrate_logit_diff / 2 / slope_agingrate,
         .keep = "unused") |> 
  summarise(across(axis,
                   list(mean = mean,
                        sd = sd)),
            axis |> 
              quantile_agingrate() |> 
              rename_with(\(x) str_c("axis_", x)),
            .by = city_code) |> 
  mutate(group = case_when(between(axis_mean, -Inf, 2035) ~ "--2035",
                           between(axis_mean, 2035, 2055) ~ "2035--2055",
                           between(axis_mean, 2055, Inf) ~ "2055--") |> 
           factor(c("--2035", "2035--2055", "2055--")))
vertex_agingrate <- agingrate_logit |> 
  left_join(axis_agingrate |> 
              select(city_code, axis_mean),
            by = join_by(city_code)) |> 
  mutate(vertex = agingrate_logit - slope_agingrate * (year - axis_mean) ^ 2,
         .keep = "unused") |> 
  summarise(across(vertex,
                   list(mean = mean,
                        sd = sd)),
            vertex |> 
              quantile_agingrate() |> 
              rename_with(\(x) str_c("vertex_", x)),
            .by = city_code)
高齢化率ピーク年のヒストグラム
axis_agingrate |> 
  ggplot(aes(axis_mean, after_stat(density))) + 
  geom_histogram(binwidth = 1) +
  geom_vline(xintercept = mean(axis_agingrate$axis_mean),
             color = "red",
             linetype = "dashed") +
  scale_x_continuous("高齢化率がピークとなる年",
                     breaks = seq(2020, 2120, 10)) +
  scale_y_continuous("密度")

高齢化率ピーク年のヒストグラム

高齢化率ピーク年別の高齢化率の推移
line_agingrate <- axis_agingrate |> 
  left_join(vertex_agingrate,
            by = join_by(city_code)) |> 
  summarise(across(c(axis_mean, vertex_mean),
                   mean),
            .by = group) |> 
  expand_grid(year = seq(2010, 2080, 1)) |> 
  mutate(agingrate = gtools::inv.logit(slope_agingrate * (year - axis_mean) ^ 2 + vertex_mean)) |> 
  select(!c(axis_mean, vertex_mean))

agingrate |> 
  left_join(axis_agingrate |> 
              select(city_code, group),
            by = join_by(city_code)) |> 
  ggplot(aes(year, agingrate)) +
  geom_boxplot(aes(group = cut_width(year, 5))) +
  geom_line(data = line_agingrate,
            color = "red",
            linetype = "dashed") +
  scale_x_continuous("年",
                     limits = c(2010, 2080)) +
  scale_y_continuous("高齢化率",
                     labels = scales::label_percent()) +
  facet_wrap(~ group,
             nrow = 1L)

高齢化率ピーク年別の高齢化率の推移

コード
kable_agingrate <- function(data) {
  knitr::kable(data, "html",
               table.attr = 'border="1"') |> 
    HTML()
}

label_axis_agingrate <- axis_agingrate |> 
  select(!c(axis_sd, group)) |> 
  mutate(across(c(axis_mean, axis_lower, axis_upper),
                round)) |> 
  rename(`ピーク年(平均)` = axis_mean,
         `ピーク年(下限)` = axis_lower,
         `ピーク年(上限)` = axis_upper) |> 
  nest(.by = city_code,
       .key = "label_axis") |> 
  mutate(label_axis = label_axis |> 
           map(\(data) {
             tags$div(
               tags$h6("高齢化率がピークとなる年(簡易推定)"),
               kable_agingrate(data)
             ) |> 
               as.character() |> 
               HTML()
           },
           .progress = TRUE))

label_vertex_agingrate <- vertex_agingrate |> 
  select(!vertex_sd) |> 
  mutate(across(c(vertex_mean, vertex_lower, vertex_upper),
                \(x) {
                  x |> 
                    gtools::inv.logit() |> 
                    scales::label_percent(accuracy = 0.1)()
                })) |> 
  rename(`ピーク値(平均)` = vertex_mean,
         `ピーク値(下限)` = vertex_lower,
         `ピーク値(上限)` = vertex_upper) |> 
  nest(.by = city_code,
       .key = "label_vertex") |> 
  mutate(label_vertex = label_vertex |> 
           map(\(data) {
             tags$div(
               tags$h6("高齢化率のピーク値(簡易推定)"),
               kable_agingrate(data)
             ) |> 
               as.character() |> 
               HTML()
           },
           .progress = TRUE))

label_agingrate <- agingrate |> 
  mutate(year = str_c(year, "年"),
         agingrate = scales::label_percent(accuracy = 0.1)(agingrate)) |> 
  pivot_wider(names_from = year,
              values_from = agingrate) |> 
  nest(.by = city_code,
       .key = "label_agingrate") |> 
  mutate(label_agingrate = label_agingrate |> 
           map(\(data) {
             tags$div(
               tags$h6("高齢化率の推移(社人研2018年推計)"),
               kable_agingrate(data)
             ) |> 
               as.character() |> 
               HTML()
           },
           .progress = TRUE))

admin_boundary_agingrate <- jpadminbdry::admin_boundary(2015) |> 
  sf::st_transform(4326) |> 
  mutate(city = city_code |> 
           parse_city(when = "2015-10-01") |> 
           city_desig_merge()) |>
  group_by(city) |> 
  summarise(do_union = FALSE) |> 
  mutate(label_city = str_glue("{pref_name(city)} {city_name(city)}") |> 
           map(\(x) {
             tags$h5(x) |> 
               as.character() |> 
               HTML()
           },
           .progress = TRUE),
         city_code = city_code(city),
         .keep = "unused") |> 
  inner_join(axis_agingrate |> 
               select(city_code, axis_mean),
             by = join_by(city_code)) |> 
  inner_join(vertex_agingrate |> 
               select(city_code, vertex_mean),
             by = join_by(city_code)) |> 
  left_join(label_axis_agingrate,
            by = join_by(city_code)) |> 
  left_join(label_vertex_agingrate,
            by = join_by(city_code)) |> 
  left_join(label_agingrate,
            by = join_by(city_code)) |> 
  mutate(vertex_mean = gtools::inv.logit(vertex_mean),
         label = list(label_city, label_axis, label_vertex, label_agingrate) |> 
           pmap(\(label_city, label_axis, label_vertex, label_agingrate) {
             tags$div(
               label_city,
               label_axis,
               label_vertex,
               label_agingrate
             ) |> 
               as.character() |> 
               HTML()
           },
           .progress = TRUE),
         .keep = "unused")
コード
set_view_agingrate <- partial(setView,
                              lng = 139,
                              lat = 38,
                              zoom = 5L)

pal <- colorNumeric(RColorBrewer::brewer.pal(11, "Spectral"),
                    domain = admin_boundary_agingrate$axis_mean,
                    reverse = TRUE)
leaflet(admin_boundary_agingrate) |>
  addTiles() |> 
  addPolygons(color = "dimgray",
              weight = 0.5,
              fillColor = ~pal(axis_mean),
              fillOpacity = 0.75,
              label = ~list_c(label)) |> 
  addLegend(pal = pal,
            values = ~axis_mean,
            title = "高齢化率のピーク年(平均)") |> 
  set_view_agingrate()
コード
pal <- colorNumeric(RColorBrewer::brewer.pal(11, "Spectral"),
                    domain = admin_boundary_agingrate$vertex_mean,
                    reverse = TRUE)
leaflet(admin_boundary_agingrate) |>
  addTiles() |> 
  addPolygons(color = "dimgray",
              weight = 0.5,
              fillColor = ~pal(vertex_mean),
              fillOpacity = 0.75,
              label = ~list_c(label)) |> 
  addLegend(pal = pal,
            values = ~vertex_mean,
            title = "高齢化率のピーク値(平均)") |> 
  set_view_agingrate()